Polynomial-Time Simulation of Pairing Models on a Quantum Computer 
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We propose a polynomial-time algorithm for simulation of the class of pairing Hamiltonians, e.g., 
the BCS Hamiltonian, on an NMR quantum computer. The algorithm adiabatically finds the low- 
lying spectrum in the vicinity of the gap between ground and first excited states, and provides a 
test of the applicability of the BCS Hamiltonian to mesoscopic superconducting systems, such as 
ultra-small metallic grains. 



PACS numbers: 03.67.Lx,74.20.Fg 



(N 
O 
O 
(N 

S3 



> 
O 



oo : 
o . 

o ■ 

^ : 

^— > ■ 
G ■ 

s : 
cr 



The potential of quantum computers (QCs) to pro- 
vide exponential speed-up in the simulation of quantum 
physics problems was originally conjectured by Feynman 
jl| , confirmed by Lloyd j^] , and later studied theoretically 
by a number of authors, e.g., [§ |, § § §• NMR-QC ex- 
periments performing quantum physics simulations were 
reported in [^| . Current QC technology is limited to fewer 
than 10 qubits and the testing of simple algorithms 
QCs of the next generation, with 10-100 qubits, have the 
potential to solve hard problems in quantum many-body 
theory. We show here how this observation can be ap- 
plied to the problem of simulating the class of pairing 
Hamiltonians with general, i.e., arbitrary long-range in- 
teractions. The pairing Hamiltonians are of wide interest 
in condensed matter and nuclear physics |p~C|] . An impor- 
tant example of a pairing Hamiltonian is the BCS model 
of low-T c superconductivity. We provide an algorithm for 
testing the validity of the general BCS Hamiltonians of 
finite particle-number systems, pertinent to nuclear sys- 
tems and mesoscopic condensed-phase systems, such as 
ultra-small metallic grains [O, [l^, [l4| . These grains 
provide a fertile testing ground for the BCS ansatz for 
the ground state wave function. The BCS wave function 
is a superposition of different Fermion numbers and is 
expected to be exact in the thermodynamic limit [fj"5[ . 
In contrast, in ultra-small metallic grains the number 
of states N within the Debye frequency cutoff from the 
Fermi energy is only ~ 100. A similar estimate holds 
for the number of states within a few major shells for 
medium or heavy nuclei. In systems with finite parti- 
cle number the BCS ansatz is doubtful, and at the same 
time exact numerical diagonalization of the general BCS 
Hamiltonian is impractical beyond a few tens of elec- 
tron pairs ]T^] . Various approximations have been pro- 
posed pH)] , but it would clearly be desirable to have an 
exact numerical solution for the problem. In H |6| ef- 
ficient QC algorithms were presented for simulating a 
many-body fermionic system. While the BCS Hamilto- 
nian describes a system of interacting fermions, it does 
so at the level of an effective field theory. This can be 
expressed in terms of an interacting spin system ]l5[ | , or 
parafermions ]17[ . Therefore the fermionic simulation al- 



gorithms H arc not directly applicable. Further, while a 
number of authors have recently considered simulation of 
one Hamiltonian in terms of another H , the connection of 
these phenomenological Hamiltonians to those of many- 
body condensed matter and nuclear physics is not a pri- 
ori clear. Here we clarify the correspondence by propos- 
ing an explicit and numerically exact diagonalization al- 
gorithm that is suitable for general pairing Hamiltoni- 
ans, and is directly implementable in NMR- type quantum 
computers ifLsf . More generally, with minor modifications 
our algorithm is applicable to all QCs with short-range 
exchange-type interactions, such as quantum dots [fl9|| . 
Using an adiabatic procedure, we show how to obtain 
only the low- lying energy spectrum, e.g., in the vicinity 
of the superconducting gap, with an algorithm that takes 
~ N A , instead of exponential, computational steps. The 
number of qubits we require equals the effective number 
of states N, so that a QC with ~ 100 qubits (neglecting 
overhead due to error correction) could solve a problem 
that is well out of the reach of current classical comput- 
ers. 

Mapping of Bosons and Fermions to Qubits. — Pair- 
ing Hamiltonians are typically expressed in terms of 
fermionic or bosonic creation (annihilation) operators, 
c$n ( c ™) an d Hn (bm), respectively, where \m\ — 
1,2,... , N denotes all relevant quantum numbers. E.g., 
the general BCS pairing Hamiltonian has the form: 
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c± m c± m is the number operator, and the 



matrix elements VJ, 



-m\ V \l, —I) (we impose no 



restriction on m, I) are real and can be calculated, e.g., 
for superconductors, in terms of the Coulomb force and 
the electron-phonon interaction JlO[] . Pairs of fermions 
are labeled by the quantum numbers m and — m, accord- 
ing to the Cooper pair situation where paired electrons 
have equal energies but opposite momenta and spins: 
rn = (p, t) and — m = (— p, i). These are degener- 
ate, time-reversed partners whose energies are considered 
phenomenological parameters [jl6| . The same idea is ap- 
plicable to nuclei, where effective pairings occur between 
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nucleons in time-reversed partners pO). N is an effec- 
tive state number, which equals the number of qubits in 
the algorithm below. E.g., in the case of metallic grains 
N is twice the the Debye frequency in units of the av- 
erage level spacing (inversely proportional to volume of 
the grain). For nuclear pairing models, N could be the 
number of states in one or more major energy shells. 

To make a connection to quantum algorithms we map 
the fermionic or bosonic operators to qubit operators. 
We denote the raising and lowering operators for the 
m th qubit by the Pauli matrices <r^, acting non-trivially 
only on the m th qubit. A "number operator" is n m = 
(af n + l)/2, where n m = 1 (0) if the m th qubit is in state 
|1) (|0)); n — ^2 m n m is the number of l's in a compu- 
tational basis state (a ket of a single bit-string), and will 
correspond, e.g., to the number of Cooper pairs in our ap- 
plications below. The computational ground state |0) = 
IO1O2 ■ • - Oat) acts as a vacuum state: n |0) = <r~ |0) = 0. 
Now we can consider three generic pairing cases and map 
them to qubits. In each case we identify fermionic or 
bosonic operator pairs that satisfy the commutation rules 
of si (2) = {cr+ , cr~ , ct^} (see |l7j for details). These cases 
are: (i) Fermionic particle-particle pairs (e.g., Cooper 
pairs): sl(2) = {c_ m c m , 4nCL m , nf n + n^ m - 1}, pro- 
vided = n_ m (a condition satisfied by Hbcs), an d 
|0) = \0) F . (ii) Fermionic particle-hole pairs (e.g., 
excitons): sl(2) = {c*_ m Cm, cj„c_ m , - n^ m } , provided 

n ra + nF m = 1 &nd |0) = J_ N • • • d 2 cLi 1 0} F . (Hi) 

Bosonic 'particle-hole' pairs (e.g., dual-rail photons in 
the optical quantum computer proposal [p0|): sl(2) = 
{b f _ m b m ,bl n b- m ,n B -n B m }, provided n^+n B m = 1 and 
|0) = b_ N - ■ ■b_ 2 b_ 1 \0) B . The three conditions above 
each restrict the dynamics to a different subspace of the 
entire Hilbert space. The conditions play the role of 
conserved quantities and only Hamiltonians that satisfy 
them preserve such subspaces. 

It is now clear how to express Hbcs in terms of qubit 
operators. In fact, a more general Hamiltonian, that is 
applicable to all cases (i)-(iii) is: 

N N 
m=l r=±Z>m=l 

(1) 

where e m = e m + V+ m and V ml = for H B cs; l,m 
now denote both state indices and qubit indices. Fur- 
ther, in the BCS case the qubit state space Hp — 
Span{|0) , ct+ |0) , <7j ; <t+ |0) , ■ • • } is mapped into a sub- 
space of the total fermionic Hilbert space where = 
n^ m . Hbcs conserves the total number operator n (the 
number of Cooper pairs). In terms of qubits, this means 
that the number of 1 1) 's in a general iV-qubit state is fixed 
by Hbcs- Thus the Hilbert space splits into invariant 
subspaces with dimension ( ) for fixed n. The prob- 
lem is reduced to diagonalizing separate blocks of size 



(^). For half-filled states in a system with TV = 100, an 
exact solution could require diagonalizing a 10 29 x 10 29 - 
dimensional matrix. Such a task is clearly unfeasible on 
a classical computer. 

Simulation of H p . — For concreteness and direct con- 
tact with feasible experiments, we limit our discus- 
sion of the simulation of H p to the nearest-neighbor 
Ising-type Hamiltonian of NMR: i?NMR = + 
X^i* Jl a f a f+ii supplemented with external single qubit 
operations F = J2u=i fi a f + fi a i- The same Hamilto- 
nian describes, e.g., a QC implementation using coupled 
Josephson junctions We emphasize that this simu- 
lation is also directly implementable in systems that use 
exchange-type interactions, since the logical operations 
for those systems are equivalent (up to polynomial over- 
head) to those using the Ising coupling |l7| . We shall 
for simplicity only explicitly discuss the case = 0, but 
the same procedure will apply also to the case of V~ l{ ^ 
(since the two cases are related by a simple unitary trans- 
formation) . From now on we denote = V m i . 

Below, we develop an explicit polynomial-time algo- 
rithm for simulating {?7 p (fcr) = exp(— iHpkr)} k !_i (t, 
T are defined later). This sequence can be Fourier- 
transformed and the spectrum of H p found [Q. How- 
ever, although this may be achieved directly using NMR 
methods, we are primarily interested in the low-lying 
spectrum (e.g., in the BCS case, near the superconduct- 
ing gap). Our algorithm therefore includes an adiabatic 
component, that allows us to probe just this part of the 
spectrum. Let us now outline the main steps in our algo- 
rithm for simulating H p using i?NMR and F. (i) Prepare 
a computational basis state \x n ) with fixed n (number 
of |l)'s). This step is well-known and needs no further 
explanation Jig] , (ii) Quasi-adiabatically evolve \x n ) to 
|'0(O))o = |<7n) + $l e n) : an approximate ground state of 
Hp (\g n ) is an exact ground state, \e n ) is a first excited 
state and 8 -C 1), with the same n as \x n ). (iii) Rotate 
|^(0)) to IV'(O)) = \g n ,n±i) + 8'\e n ,n±i), a state that in- 
cludes contributions from n ± 1 as well, (iv) Implement 
U p (t) = exp(-iHpt) on \tp(0)). (v) Measure. Repeat 
steps (i)-(v) while increasing t in step (v). We describe 
each of these steps in detail, starting for simplicity from 
step (iv). 

Step (iv): Implementation of exp(—iH p t). — In NMR 
one can only control ff (or ff) directly, while all ui, Ji 
are always on |l8| . Also, Ji usually is positive. A power- 
ful method that allows us to deal with such constraints 
(that are not unique to NMR) is recoupling (e.g., p2| ). 
The idea is based on elementary angular momentum the- 
ory. We define C\ o e l6B = e ifA e ieB e -i 9 A^ where A B 
are generators of su(2) (e.g., two Pauli matrices), and/or 
{A, B} = while A 2 — 1. This recoupling sequence can 
be interpreted as the application of time-reversed pulses 
(e ±lipA ) before and after periods of free evolution e l0B . 
Special cases of interest are (i) C^ 2 o e l6B = e~ l8B , 
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ing 0(N 3 ) steps. However, H p also contains the term 

Hq =J2iLi % a i> wmc h d° es n °t commute with Hi. 
Clearly, by turning on single qubit NMR erf terms for 
times T\ so that loiti = bit, we can simulate Hq directly 
using N steps. The non-commutativity implies that we 
need a short-time approximation in order to simulate the 



-iHf full U p (t) = exp(- 
U P (r) 



iH p r) 



—iHoT—iHiT 



+ 0(r 2 ). 



(2) 



FIG. 1: Quantum circuits to simulate e~ lH ° r (a) and e~ lH ' T 
(b) for the two qubit case. Time flows from left to right. 
X = a x . The recoupling procedure yielding U z (t) — 
exp(— iJ\Tv2.a\ erf) is in the box in (b), and is repeated with- 
out detail. We set uiiTi = £;t (i = f , 2) and U\t\2 = | V12 1 t. 
Rectangular boxes connecting two qubits denote evolution un- 
der I/nmr for the indicated time. 

(ii) C^ /4 o e l9B = e id ^ BA l Thus, to obtain evolu- 
tion under ^crf we apply the (unoptimized) recoupling 
sequence exp(-if-aft) = ( e -iH NUR t/4 Tie -iH NUH t/4 T /y ^ 
where T; = ®j^;crj, T[ = (gi^c? where the prime indi- 
cates that j is even (odd) if I is even (odd). This takes 
3N pulses. Fig. 1(a) illustrates an optimized circuit for 
N = 2. Similarly, we can evolve under any term cr|cr| +1 
using ~ 77V recoupling steps. 

Next, we need to show how to simulate long-range 
interactions using -Hnmr and F. The set {Xi m = 

\ Of^m + °f °m) , Y lm = \ (of '<?%, ~ <jf (J y m ) , Z \ m = 

i (erf — °~m)} forms an su(2) algebra, and commutes with 
a^ + af for any I, m Thus Cj/ 2 +i o Z u+1 = -Z tt i + i, 
while C^ 2 +i o (erf + of +1 ) = (of + crf +1 ). Adding yields 
<%£ + A°U°t) =af_ 1 af +1 , so that C^oe**?-^ = 
e *$<r*_ 1 of +1 xhus CZ/ 2 acts as a nearest-neighbor ex- 

change operator. In order to implement C X[[+i using 
-Hnmr and F note that: 

a i +cr i+i CT i ; + CT !+i 

It is simple to check that to create all possible couplings 
crfer^ in this manner requires 0(N 3 ) steps. This pro- 
cedure allows us to use the short-range NMR Hamil- 
tonian to simulate Jiafa^ with \l — m\ arbitrary. Let 
us now show how to turn this into a simulation of 

H I = I Eilm=l V ml(Vm a f + °"m°f )• Suppose that H p 

evolves for time r. We can turn on —Jiafa^ for a 
time T m i such that 2Jit„u — \V m i\T (for a BCS Hamil- 
tonian V m i < 0). Doing this for all couplings sep- 
arately (in series) shows that the evolution operator 
U z (t) = exp(-| Y^i> m v im <J i <J m T ) is obtained using the 
same 0(N 3 ) steps. By adjusting single-qubit opera- 
tion times, we can implement U a = exp(i j J2i of), to 
yield: exp(-i-ffjT) = (U'^U z {t)U x ) (U v U z {t)U v ^), us- 



When the additional recoupling steps needed to turn 
off unwanted interactions (which we ignored above) are 
taken into account, using the method of p2| , we find 
that C/ p (r) requires a total of s(N) = -§A^ + 22. N - 
47 ^y3 _|_ 28 s t e p S> This result may be improved some- 
what if parallel operations are allowed. E.g., in Fig. 1 
we show optimized circuits implementing eT ° T and 
e~ IT for N = 2 qubits. If Hnmr contains beyond- 
nearest- neighbor interactions then at most 0(N 5 ) steps 
are needed. The effect of the 0(t 2 ) errors in quantum 
algorithms due to the short-time approximation has been 
analyzed, e.g., in H. By concatenating short-time evolu- 
tion segments one can then obtain the finite time (fcr = i) 
evolution operator U p {t) » (U p (T)) k j|, in a total of 
ks(N) steps. 

Step (ii): Adiabatic Evolution. — Let 2A be the gap be- 
tween the ground and the first excited states, and let < 
c(t) < 1, c(0) = 0, c(T) = 1, be a slowly varying function, 
i.e., 2n/T < 2 A (e.g., c(t) = t/T). Consider the time- 
ordered evolution U a d{t) — Texp(— i J Q * H(s)ds) under a 
time-dependent Hamiltonian H(t) = Hq + c(t)Hj. For 
sufficiently small r this factors into a product 

C/ad(fcr) » e - iH ^ T ■ ■ ■ e -iH{2r)r e -iH{r)r + 

(3) 

where exp(— iH(Jt)t) rs exp(— iHqt) exp(— ic(jr)HiT) 
(j = 1, k), and now we choose times T m i(j) (for turn- 
ing on -Jiafa^) such that 2J ; r mi () = \V m i[r c{jr). 
Since c{t) is slow, U a d(kr) will represent an adiabatic 
evolution. The adiabatic theorem then ensures that the 
system will be in an eigenstate of H p = H(T) at T = kr, 
provided the initial state is in an eigenstate of H . More- 
over, this will be a ground state \g n ) of H p (a state with 
fixed n) if the initial state is the ground state of Hq (a 
computational basis state \x n )) 0]. In order to probe 
the low-lying spectrum we may slightly relax the adia- 
batic condition ir/T < A, or S; > tt/(tA). This can 
be defined in terms of the adiabatic expansion where the 
first order constraint is the usual adiabatic assumption. 
Here we only wish to satisfy the second order condition 
p5[ . Then we obtain a state |i/ , (0))o ~ \g n ) +^|e n ) which 
contains a small (6 -C 1) component |e n ) of some of the 
low-lying excited states of H p (with the same n). 

Steps (iii),(v): Measuring the Spectrum. — In NMR one 
measures the free-induction-decay (FID) signal, given 
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by V a (t) oc Tr(p(t)a~), where p(t) is the system den- 
sity matrix and a is the index of the measured spin 
(qubit) pl|. To probe states with different n, we ro- 
tate to |^(0)) = e- fawr S|^(0)) PS |ff«,«±i> + 0'\e n ,n±i), 
where 0',lu <C 1, a state that includes contributions 
from n ± 1 as well [siep fni/]. This is simple to do 
using the method of step (iv). Combining steps (ii)- 
(iv), we have p(t) = U p (t)\ip(0))(ip(0)\Ul(t). To re- 
late V a (t) to the spectrum of the pairing Hamiltonian 
we introduce an appropriate basis. A complete set of 
conserved quantum numbers are the number of Cooper 
pairs n (= the number of l's in a computational ba- 
sis state, lowered by cr~), the energy E n> i for fixed 
n, and a state degeneracy index fy. Thus our basis 
states are labeled by \n,i,/3i) and p(t) can be expanded 



as J2 B n.i,p t B, 



i(E„ 



-B„,i)t, 



with 



Pi). We have 



V a (t) <x 



m.n ij 



m.j;n,i 



(4) 



where C 



(a) 

m,j;n,z 



^(mJiPjWa \n,i,Pi) 
Fourier transforming, we obtain the energy 



ft/3. 



m,3 



spectrum S(u) = J2n,i,j Cn-l,j;n,i 5 { w ~ ( E n-i,j - En,i)), 
with the gap defined as 2A n = E Ui i — E n< o. Ideally, 
A„ can be found from a few runs with different initial 
n. There are two complications in practice: (i) Finding 
A n in this manner depends on the coefficients C^lj -. n i 
not vanishing. By measuring all qubits a, it is likely 
that sufficiently many non-zero coefficients will be 
available, (ii) The sharpness of the 8 functions depends 
on how densely the signal V a (t) is sampled. To resolve 
the gap, we will need to sample with a resolution 
Aui — 2ir/T < A„. Recall that -ffecs conserves n. 
Thus the number of r-intervals required for fixed n is 
k(n) 3> 7r/(rA n ), which is just the adiabatic condition 
again. A total of ifc(n) 2 elementary evolutions steps, 
each simulating evolution under H p for length r, will 

T It 

thus be needed to simulate {U p (kT)} k L. 1 , and each 
such step takes s(N) logic gates. The longest single 
run takes k(n)s(N) steps, while ^k(n) 2 s(N) is the 
total run-time of the algorithm, if the algorithm is to 
succeed in the absence of error correction, then we must 
have k(n)s(N) < T 2 /ri ogic , the ratio of decoherence 
to logic gate time. For NMR, T 2 /ri og i c can be ~ 10 5 . 
To estimate k(n) we need r and A„. The gap can be 
estimated experimentally, for nuclear and BCS systems 
using material dependent parameters |l^, 0|. Recall 
that r is related to the short-time approximation which 
allowed us to neglect commutator terms in the expansion 
of U ad (t). Since e^ A+B ^ » e Ar e Br e ~^ A ^ T \ we need 
to estimate when \[A, B]t\ -C mm(|A|, \B\). To obtain 
a rough estimate we consider a reduced BCS model 
|L4j: Vmi = -V < 0, ei = £ + Id. In the BCS case the 
level spacing d C V, but e ^ V. Letting A = eiaf, 
B = VX lm , we have \[A,B]\ = \V (ej - e m ) Y lm \ > Vd, 



while min(|^4|, \B\) = V. Thus the short-time approxi- 
mation is valid when r <C \jd. Using k(n) ^> 7t/(tA„) 
and s(N) w 9N 4 we thus have k{n) s(N) > 30^- A^ 4 . In 
the BCS case d/A n <C 1. Assuming d/A n — 0.1 we find 
k(n) s(10) > 3 x 10 4 , so that a simulation with N < 10 
qubits seems to be within the reach of present day NMR 
simulations pj| . 

In order to illustrate the algorithm, consider a 
simple example, the circuit for which is given in 
Fig. 1. When N — 2 the computational basis states 
are: {|00) , |01) , 1 10) , |11)}, with n = 0,1,1,2 Cooper 
pairs, respectively. Diagonalizing H p yields the en- 
erg y spectr um: {E n } = {^o = -(ei + e 2 )/2,Ef = 
±Ve 2 + V 2 ,E 2 = (ei +e 2 )/2}. Steps (ii)-(v) of the algo- 
rithm can be carried out analytically. Fourier transform- 
ing the FID signal yields four spectral lines from which, 
e.g., the n = 1 gap can be found as 2Ai = E^ — Ey . 

Conclusions. — We have proposed an efficient algo- 
rithm for finding the low-lying spectrum of pairing mod- 
els with arbitrary long-range interactions, such as the 
BCS Hamiltonian. This establishes a link between quan- 
tum computers (QCs) of the next generation (10-100 
qubits) and outstanding problems in finite-system quan- 
tum physics, such as the applicability of the BCS model 
to mesoscopic solid-state and nuclear systems. It would 
be interesting to implement the algorithm using current 
NMR-QC know-how, thus extending the experimental 
repertoire of QC physics simulations ||. 
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